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ABSTRACT 

Wider use of ceramic matrix composites (CMC) will require the development of advanced structural 
analysis technologies. This report focuses on the use of an interactive model to predict the time-independent 
reliability of a component subjected to multiaxial loads. The deterministic, three-parameter Willam-Warnke 
failure criterion serves as the theoretical basis for the reliability model. The strength parameters defining the 
model are assumed to be random variables, thereby transforming the deterministic failure criterion into a proba- 
bilistic criterion. The ability of the model to account for multiaxial stress states with the same unified theory is 
an improvement over existing models. The new model has been coupled with a public-domain finite element 
program through an integrated design program. This allows a design engineer to predict the probability of 
failure of a component. A simple structural problem is analyzed using the new model, and the results are com- 
pared to existing models. 


CHAPTER I 
INTRODUCTION 

The ability of structural components, fabricated from both monolithic and composite ceramic material 
systems, to maintain their structural integrity while subjected to thermomechanical loads is beginning to capture 
the attention of many design engineers. Attractive properties such as low density, high strength, high stiffness, 
creep resistance, and corrosion resistance are allowing ceramic materials to supplant metal alloys in numerous 
applications. Current applications include heat exchangers, cutting tools, and wear parts. Larsen and Vyas 
(1988), Buljan, Paslo, and Kim (1989), and Clarke (1990) present commercial data regarding the expanding use 
of ceramic components for these rigorous applications. Unlike some metal alloys used in demanding service 
conditions (notably the superalloys), ceramic components are fabricated from nonstrategic materials. This has 
helped spur research efforts in both processing technology and structural analysis. This report focuses on issues 
related to the field of structural analysis, where design protocols are replacing the ad hoc trial-and-error method 
of developing and testing structural prototypes. 

In the field of material science, efforts to improve the structural performance of ceramic materials include 
adding a second ceramic phase to the matrix. This second phase can take the form of whiskers, short (usually 
chopped) fibers, continuous fiber reinforcement, and woven fabrics. This report will focus on ceramic com- 
posites that incorporate whisker reinforcement (and under certain conditions, particulate reinforcement). The 
addition of whiskers improves the failure behavior of the material system by arresting crack growth in the 
matrix by pinning, bridging, and deflecting cracks. The improvement of fracture toughness, usually in certain 
material directions, is dependent on processing. As a result, this material can exhibit anisotropic behavior. How- 
ever, if the whiskers are homogeneously distributed and randomly oriented, the isotropic nature of the matrix 
material is preserved. The work presented here will deal exclusively with the isotropic whisker-toughened 
material system. Analytical efforts that allow for material anisotropy are mentioned in chapter V. 

Even though the second phase enhances the failure behavior of the material, whisker-toughened ceramics 
still fail in a brittle fashion. In addition, there is a great deal of intrinsic variability in the strength of this 
material. Failure of structural components fabricated from whisker-toughened ceramics is governed by random 
flaw populations inherent to the material’s microstructure. Usually these material imperfections are generated 
during processing. It is assumed that the location and orientation of the flaws are randomly distributed 



throughout a component. The resulting scatter in failure strength of these materials requires a departure from 
traditional design philosophies. The random nature of the microstructural flaws forces the design engineer to 
rethink the design philosophy that treats material strength as a single-valued design parameter. For monolithic 
ceramics, the factor of safety approach (a deterministic design procedure commonly used for metal alloys) has 
been abandoned in favor of a reliability-based approach. Work by Gyekenyesi (1986), Cooper, Margetson, and 
Humble (1986), and Lamon (1990) are representative of the reliability design philosophy used in analyzing 
structural components fabricated from monolithic ceramics. 

Adopting a similar probabilistic philosophy for the structural analysis of a component fabricated from 
whisker-toughened material allows the design engineer to account for brittle behavior, variability in strength, 
and decreasing bulk strength with increasing component volume (the so-called size effect). Using probabilistic 
methods, the component is discretized using finite element techniques, and each discrete element is treated as a 
link in a chain. Philosophically, this means that when one element fails, the component fails. Thus the com- 
ponent is only as strong as the weakest link in the chain. From the standpoint of reliability theory, the compo- 
nent is treated as a series system, where failure of the system occurs when one of the subsystems fails. 
Alternatively, in a parallel system, failure of a single subsystem does not cause the system to fail since the 
remaining elements of the system may sustain load through redistribution. Models that use the analogy of a 
parallel system lead to what has been referred to in the literature as bundle theories. The basic principles under- 
lying bundle theories were originally discussed by Daniels (1945) and Coleman (1958). Bundle theories have 
been applied exclusively to long-fiber ceramic composites. Thus, further discussion of these theories will not be 
pursued here. See the work of Harlow and Phoenix (1981), and Phoenix (1974, 1979) for an in-depth treatment 
of the bundle theory. 

In general, two categories of weakest link theories have emerged. One group is based on the principles of 
fracture mechanics. The other group adopts a phenomenological viewpoint. The fracture mechanics approach 
assumes that the stress state in the near vicinity of the critical crack and the orientation of the crack are the con- 
trolling design variables. Material strength and crack orientation are treated as random variables. All other 
design variables (e.g., load, geometry, stiffness, etc.) are treated in a deterministic fashion. In contrast, phenom- 
enological reliability models take a more global approach. Only material strength is treated as a random variable 
since attention is not focused on a critical flaw. Phenomenological models can be either interactive or nomnter- 
active. Interactive models allow functional forms that include terms that are products of different material 
strengths. Noninteractive theories allow material strength parameters to appear only as separate and distinct 
terms. Throughout this report, the fracture mechanics and the phenomenological criterion will be discussed for 
the purpose of comparison, but the attention will be focused mainly on interactive models. In chapter II, a litera- 
ture survey is presented that outlines different methods of modeling reliability. 

With the exception of the work by Adams and Sines (1978), Alpa (1984), and Powers (1989), reliability 
theories for ceramic components have neglected compressive stress states and the effect of hydrostatic stress in 
particular. Models such as the principle of independent action (PIA), which was originally proposed by Barnett 
et al. (1967) and Freudenthal (1968), and the familiar Batdorf theory (Batdorf and Crose, 1974) do not allow 
compressive states of stress to influence component reliability. Since the compressive strength of ceramic 
materials is often an order of magnitude larger than the tensile strength, compressive stress states were assumed 
not to contribute to failure, or treated in an ad hoc fashion in a manner similar to Gyekenyesi (1986). Although 
data in the open literature are limited, experimental evidence by Adams (1975), Ikeda and Igaki (1984) and 
Ikeda. Igaki, and Kuroda (1986) clearly indicates that compressive stress states have a decided effect on ceramic 
materials. The phenomenological criterion that is discussed later allows for multiaxial states of stress, and 
specifically treats compressive stress in a rational manner as outlined in chapter IV. The criterion (and the 
reliability model that is derived from this criterion) is unified in the sense that ad hoc rules are not used to 
model different regions in the stress state. The analytical details of the parent deterministic failure model are 
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given in chapter III. Failure surfaces projected into various stress spaces are presented to illustrate different 
mechanistic aspects of the theory. Finally, the tests necessary to determine the parameters are outlined. 

Casting the deterministic failure theory into a reliability model using Monte Carlo methods is presented in 
chapter IV. Numerical aspects of the Monte Carlo simulation are discussed. Features of the interactive reliability 
model are compared with existing models. The reliability model is incorporated into a test-bed software program 
given the acronym TCARES (Toughened Ceramics Analysis and Reliability Evaluation of Structures), which 
was originally discussed by Duffy et al. (1989). Coupling the reliability algorithm with a general-purpose finite 
element program (i.e., MSC/NASTRAN) enables one to predict the time-independent reliability of a structural 
component. A structural component is analyzed that illustrates the interactive model highlighted in this report. 
These results are compared with an analysis made using previous models that did not allow compressive stress 
states to affect component reliability. It is shown that these previous models yield unconservative results in cer- 
tain situations. 

Chapter V summarizes this effort and indicates future research. Future direction includes improving 
numerical efficiency through the use of fast probability integration techniques proposed by Wu (1984). Applying 
this type of analysis to anisotropic whisker-toughened ceramics is outlined. 
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CHAPTER n 


SURVEY OF RELIABILITY THEORIES 


Traditional failure analyses of structural components have used deterministic approaches where failure is 
assumed to occur when some allowable stress level, or equivalent quantity, is exceeded. This assumes that 
deformation is not controlling component design. Since structural ceramics maintain high stiffness, even at 
elevated temperatures, deformation has not played a significant role in component design. Certain design 
methods have attempted to incorporate the relevant physics of failure using fracture mechanics. Here the critical 
design parameter is the stress intensity factor, which takes into account load and component geometry. In this 
approach the stress intensity factor is compared to a fracture toughness value that is a characteristic property of 
the material. However, for most structural ceramics the combination of ultimate strength and fracture toughness 
(quantified by K 1C ) yields flaw sizes so small that current nondestructive evaluation (NDE) methods are unable 
to detect the critical defect. On the other hand, phenomenological failure theories make use of macroscopic 
strength parameters that do not focus on a critical microstructural defect. Multiaxial failure theories can be sys- 
tematically formulated using this approach if the material is homogeneous, with strength properties that can be 
deduced from well chosen phenomenological experiments. Failure theories such as the maximum normal stress, 
the maximum normal strain, the maximum shear stress, and the maximum distortional energy criteria are 
examples of phenomenological models that are successful in predicting the onset of brittle failure or yielding. 
However, for reasons mentioned in the introduction, these deterministic techniques are not relevant when 
analyzing structural components fabricated from ceramic-based material systems. 

Weibull (1939, 1951) proposed the first probabilistic model that accounted for scatter in failure strength 
and the size effect encountered in brittle materials. His approach is based on the weakest link theory (WLT) 
attributed to Midgley and Pierce (1926). This earlier research (sponsored by the textile industry) focused on 
modeling yarn strength. Unlike Midgley and Pierce, who assumed a Gaussian distribution for yarn strength, 
Weibull proposed a unique probability density function for failure strength that now bears his name. Weibull’s 
two-parameter probability density function has the following form: 
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for a continuous random variable x, when x > 0, and 

f(x) - 0 


(2.1) 


( 2 . 2 ) 


for x < 0. The cumulative distribution function is given by the expression 


F(x) 


= 1 - exp 



for x > 0, and 


F(x) - 0 


(2.3) 
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for x < 0. Here a (> 0) is the Weibull modulus (or the shape parameter), and (? (> 0) is the scale parameter. 
Reliability theories with theoretical frameworks based on Weibull’s original concepts are presented in this 
chapter. Theories based on phenomenological principles and fracture mechanics theories are discussed. Initially, 
Weibull’s (1939) normal stress averaging technique is presented. This is followed by a discussion of the PIA 
model, and recent extensions of the PIA model to composite materials. Next, a reliability model developed by 
Batdorf and Crose (1974), founded on principles of fracture mechanics, is presented. Finally, a model that 
accounts for compressive states of stress (Powers, 1989) is discussed. 


Weibull’s Normal Stress Averaging Method 

Weibull adopted the weakest link theory where a brittle material is considered a chain with links connected 
in series. The overall strength of a brittle component is then governed by the strength of its weakest link. 
Focusing on a single link, the failure probability of an individual link can be expressed as 

p f - ipAV < 2 - 5 ) 


where AV is the volume of the link, and ip is a failure function per unit volume of material. By defining r 
as the reliability of a single link, then 


r = 1 - tpAV 


( 2 . 6 ) 


The failure of an individual link is assumed to be an independent statistical event, implying that the events lead 
ing to failure of an individual link are not influenced by other links in the chain. As a result, the reliability of 
the component, denoted as R, becomes 


R = lim 

1 

= lim 
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n 
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N-»°° 




where N is the number of links in the component, x^) is the failure function per unit volume at position 

Xj within the component, and a y is the Cauchy stress tensor at x^. Unless noted otherwise, lowercase Roman 
letter subscripts in italics denote tensor indices with an implied range from 1 to 3, and Greek letter subscripts 
are associated with products or summations with ranges that are explicit in each expression. By adopting the 
argument originally proposed by Weibull, the reliability of the component takes the integral form 

R = exp ["Jv M, ( <7 y’ x, ) dV ] (2 ' 8j 

where V is the volume of the component. Similarly, the probability of failure for the component takes the form 

P f = 1 - R - 1 -cxp[-J^tp(a, 7 , x ( )dvj (2 ' 9) 

The state of stress in every’ link of the chain must be characterized to conduct a reliability analysis using 
equation (2.8). This approach lends itself to analytical techniques that use finite element methods. If a 
component is modeled as a chain with individual links connected in series, then each link would correspond to 
an element within a mesh. The reliability models that are discussed in this section, and in later sections, adopt 
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this theoretical framework (eqs. (2.5) to (2.9)) in computing component reliability, and all are amenable to finite 
element methods. The differences between various reliability models occur in the formulation of the failure 
function tp. 

Weibull assumed that the failure strength of a specimen subjected to a uniaxial state of stress is a random 
variable. Application of equation (2.9) for a uniaxial tensile stress field in a homogeneous isotropic material 
yields 


Pf 



( 2 . 10 ) 


where a is the applied tensile stress, and V is the volume of the specimen. For this case Weibull took the fail- 
ure function (tp) as 


ip - 


/ '<y' vct 
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( 2 . 11 ) 


Here k is referred to in the literature as Weibull’s coefficient for a uniaxial state of stress. Weibull extended 
this uniaxial model to multiaxial states of stress by defining an average tensile stress. He defined this average 
tensile stress by considering the stress traction on an arbitrary plane (see fig. 2.1). Specifically the shear 
component of the stress traction r is ignored, and it is assumed that only the normal component o causes 
failure. To gain a clear understanding of this method, consider a sphere centered at the origin of the coordinate 
axis (Xj, x 2 , Xj) associated with the principal directions Oj, a 2 , and a 3 . The normal component of stress 
traction acting on an arbitrary plane is given by the expression 

a = sin 2 4>(aj cos 2 0 + a 2 sin 2 ©) + a 3 cos 2 <J> 

Here (j) and 0 are the polar and azimuthal angles of the unit vector normal to the arbitrary plane in stress 
space. These angles and their relationship to the coordinate axes (Xj, x 2 , x 3 ) are shown in figure 2.2. Since the 
normal stress component varies with each planar orientation defined by <f> and 0, Weibull defined the follow- 
ing weighted average 



8 j r j;Vsin 4> d0 d<J> 


(2.13) 


Here A is the area of the unit sphere, and the limits of 0 coincide with those orientations where a changes 
from a tensile stress to a compressive stress. For the limiting case where a is tensile over the entire sphere, 0 
takes the value of n/2. Both limits are zero if a is compressive over the entire unit sphere, and this results in 
(d) a = 0. The general convention analyzes this integral over one octant of the unit sphere, which accounts for 
the limits of integration for (j> (i.e., 0 and n/1), and the factor 8 in the numerator. 
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The quantity ( o) a given by equation (2.13) is equated to o a in equation (2.1 1) such that 


(2.14) 


ip = k(d) a 


where k is Weibull’s coefficient for a multiaxial stress state. This coefficient is the multiaxial extension of the 
parameter k defined in equation (2.11), and must be defined in a consistent fashion such that equation (2.14) 
yields equation (2.1 1) for a uniaxial state of stress. Equating Weibull’s multiaxial formulation of ip defined in 
equation (2.14) to the uniaxial case defined in equation (2.1 1) results in the following relationship (see 
Gyekenyesi (1986) for details): 

k - k(2« + I) - (2 ° - *> (2 '5) 

S” 


Thus for the multiaxial state of stress. 


P f = 1 -exp k(o ) a dV 


from which the uniaxial form expressed in equation (2.10) can be obtained. Although the extension to multiaxial 
states of stress described here is intuitively plausible, it is somewhat arbitrary because it disregards the shear 
component of the stress traction. In addition, since the method lacks a closed-form solution, use of this model 
requires computationally intensive numerical methods. 


Principle of Independent Action Method 

Barnett et al. (1967) and Freudenthal (1968) proposed an alternative to Weibull’s normal stress averaging 
approach for multiaxial states of stress. Here only principal stresses are considered, and the basic assumption is 
that each acts independently in reducing the survival probability of an element (hence the name principal of in- 
dependent action). The failure function for this theory takes the form 

a 3 (2.17) 

T. 



where > o 2 > o 3 > 0. Since the principal stresses appear in separate terms (i.e., they do not interact) in the 
formulation of ip, this model is classified as a noninteractive reliability model. Qualitatively, the PIA theory is 
equivalent in a probabilistic sense to the maximum stress failure theory. 


Extensions of the PIA Method 

Duffy and Arnold (1990) formulated an extension of the PIA model for transversely isotropic materials. A 
unit vector was used to identify the local material orientation, and, subsequently, to define stress invariants. The 
unit vector d ( was introduced to define the direction normal to the plane of isotropy. Here the failure function 
ip depends on the stress state and local material orientation such that 
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v - v (<v d ,) 


( 2 . 18 ) 


Since ip is a scalar valued function, it must remain form invariant. To ensure this, an integrity basis was 
developed for ip that contained certain combinations of invariants of the Cauchy stress tensor and an orienta- 
tion tensor defined as d ( dy. The invariants formed an integrity basis and were used to construct other invariants 
that correspond to specific components of the state of stress in an element. This approach yields the following 
functional dependence 


ip = vp(ti,i 2 , i 3 , U) (2I9) 

Here Ij corresponds to the magnitude of the stress vector projected onto the material orientation vector 

d-. The invariant i 2 represents the magnitude of the shear component of the stress traction. The invariants t 3 
and I 4 represent the magnitudes of the maximum and minimum principal stresses in the plane of isotropy. 

Thus each invariant corresponds to a strength in a well-defined material direction, and because of this the 
invariants can be treated as random variables with underlying Weibull distributions. Analogous to the principle 
of independent action, it is assumed that the different invariants are statistically independent such that the failure 
function takes the form 




V 

«1 

*L. 

M 


<i 3 >' 

«3 

'<l 4 >' 

Pi J 

r 

hJ 

i 


1^3 J 


[P 3 J 


( 2 . 20 ) 


where the individual a’s and P’s are the Weibull parameters associated with a strength variable in a particular 
material direction. It should be noted that compressive stresses associated with <!j>, <i 3 >, and <I 4 > are assumed 
not to contribute to a reduction in reliability such that 


<V 


ii I, > o 

o < o 


( 2 . 21 ) 


and 


<i 3 > 


4 ^3 ^3 > ^ 

0 1 3 < 0 




U U >0 

o i 4 < o 


( 2 . 22 ) 


(2.23) 


Duffy and Manderscheid (1990) formulated an extension of the PIA model for orthotropic materials. Here 
two mutually orthogonal unit vectors (a.- and b ( ) were used to define local material directions. An integrity basis 
was developed from the functional dependence 
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V = vfay.afbj 


( 2 . 24 ) 


where a,b serves as a direction tensor. Like the transversely isotropic case, the failure function depends on 
certain invariants that correspond to components of the local stress tensor, that is. 


IP 



( 2 . 25 ) 


Here the invariants I, and I 3 represent the ^magnitude of the normal stress components in the directions of a, 
and b, respectively. The invariants I 2 and I 4 represent the shear stresses across the directions a, and b, 
respectively. The invariant I 5 represents the normal stress in the direction defined by the cross product of the 
vectors a- and b. Once again, it is assumed that the invariants are statistically independent such that 
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where the individual ot's and P's are again associated with strength variables in particular material directions. 
The invariants I,, I 3 , and I 5 are normal stresses, and do not contribute to a reduction in reliability if they are 

compressive; that is. 


<Ii> 


o r t < o 


( 2 . 27 ) 



I 3 > 0 

r 3 <o 


( 2 . 28 ) 


and 



I 5 > 0 
I 5 < 0 


( 2 . 29 ) 


The reliability models for both the transversely isotropic and orthotropic materials allow the material orien 
tation to vary along a family of curves within the component. Thus the material is locally anisotropic. The. 
models were constructed using invariant formulations which indicate the maximum number and forms of the 
stress invariants necessary to define the failure function tp. In both cases, a subset of the integrity basis or ip 
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was constructed, resulting in reliability models that are similar in nature to the PIA reliability model for mono- 
lithic ceramics. 


Batdorf s Theory-Surface Flaw Analysis 

Reliability theories based on fracture mechanics assume that failure of a component emanates from a single 
flaw with a critical size and orientation. This flaw belongs to a population that in general contains surface flaws 
and volume flaws. Either type of flaw is assumed to be uniformly distributed and randomly oriented. Surface 
flaws are imperfections that are the result of machining, grinding, or other surface finishing operations. Volume 
flaws are the direct result of processing. Both ty pes of flaw populations exhibit different failure behavior charac- 
terized by distinct strength distribution parameters. For surface flaws, the presence of a traction-free surface 
reduces a three-dimensional state of stress to a state of plane stress. Because of this simplifying condition, the 
details of a surface flaw analysis are presented. In general, the surface flaw analysis can be viewed as a special 
case of volume flaw analysis. 

Batdorf and Heinisch (1978) proposed a surface flaw model which is an extension of an earlier volume 
flaw analysis proposed by Batdorf and Crose (1974). This two-dimensional theory is based on weakest link prin- 
ciples, and accommodates mode I, mode II, or mixed mode fracture criteria. For the following discussion, the 
coplanar strain energy release rate, a mixed mode criterion, is used. This criterion allows for mode I and 
mode II behavior, and takes the form 
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1C 
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II 


K 


1C 


(2.30) 


where K, and K n are the mode I and mode II stress intensity factors, respectively. These stress intensity 
factors are functions of the applied far-field stress state and the crack geometry. For a Griffith crack (a sharp 
through-crack of length 2a), the stress intensity factors for a two-dimensional infinite plate are 


Kj - Qyfita 


(2.31) 


and 



(2.32) 


Here a is the far-field normal stress, and r is the far-field shear stress. Substituting equations (2.3 1 ) 
and (2.32) into equation (2.30) yields the following expression: 
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At this point, index notation is briefly suspended in order to use notation that has been widely accepted in dis- 
cussing BatdoiTs theory. By defining 
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(2.34) 
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cr 
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as a critical stress, then equation (2.33) (together with eq. (2.34)) defines a failure envelope which is dependent 
on the material’s fracture toughness and crack size. This failure envelope is shown in figure 2.3. Since the size 
(denoted as "a" in the denominator of eq. (2.34)) and orientation of the critical crack will vary, a cr represents a 
random variable. Thus there is a family of failure envelopes corresponding to each value of the random variable 
a cr . Batdorf used this concept to transform a deterministic fracture criterion into a reliability model. 

The Batdorf theory stipulates that the probability of failure of a single link is the product of two probabili- 
ties; that is. 


Here p, is the probability that a crack exists such that the applied stress is in the range of o cr to (o cr + da cr ). 
This probability is defined by the expression 


Pi 


AA 


dN(<g 
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da„ 


(2.36) 


where AA is the differential area of the link, and N(o cr ) is a crack density function (i.e., N has units of cracks 
per unit area). The crack density function is defined as 


and quantifies the number of cracks per unit area with an applied far-field stress that is greater than or equal to 
a Here k B is the Batdorf constant (which is functionally dependent on the Weibull scale parameter) and a 
is'the Weibull modulus of the material. Both parameters can be determined experimentally by using the para- 
meter estimation techniques outlined in Pai and Gyekenyesi (1988). 


The quantity p 2 is the probability that a crack is oriented such that the critical stress is exceeded by an 
applied far-field stress. To illustrate this point, consider a state of plane stress using Mohr s circle where 
0 < a < o | (see fig. 2.3). Define the angle oo as the initial orientation of the crack where the critical stress is 
exceeded. Note that the crack orientation varies between 0 and 2 tt radians, and there may be more than one 
orientation where the critical stress is exceeded. Next, overlay the fracture envelope expressed in equation (2.33) 
on the same set of axes. The angle between the O axis and the line OA defines 2co. Here point O is the 
center of Mohr’s circle defined by the in-plane principal stresses and a 2 . Point A is the intersection of the 
failure envelope given by equation (2.33) and Mohr’s circle and defines the initial orientation of a crack where 
the critical stress is exceeded. Any point on Mohr’s circle outside the failure envelope represents a possible 
failure stress state. Thus the probability that a crack is oriented such that the critical stress is exceeded by the 
applied far-field stress is 
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where 0 <p 2 < I . 

The geometry in figure 2.3 is used to established the form of equation (2.38). The applied far-field stress is 
equal to o cr at point A; thus 


Hence, 


CO = cos 




(2.39) 


(2.40) 


when o 2 < a cr < and 


P2 " 1 


( 2 . 41 ) 


when a cr < o 2 . For the case where a 2 < < 0 cr , 


p 2 = 0 


(2.42) 


Following the argument outlined in equations (2.5) to (2.7), with the results of equations (2.36) and (2.40), 
then the reliability of a component is given as 
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Similarly, the probability of failure of a component becomes 
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Note that the limit of integration defined by A is the area of the component, and the limits of integration for 
a cr assume that Oj is the largest principal stress. In the context of the preceding discussion concerning reli- 
ability models, the failure function for Batdorf’s model is defined by the integral 
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The above analysis assumes that compressive stresses do not contribute to failure. For discussion on Batdorfs 
model for a volume analysis, see Batdorf and Crose (1974). 


Power’s Extension of Batdorf s Theory 

The Powers model (Powers, 1989) represents a hybrid approach to reliability analysis in the sense that 
Batdorfs approach is adopted to predict reliability of a single link when both of its principal stresses are tensile, 
and the Mohr-Coulomb failure theory is used when both are compressive. When one principal stress is tensile 
and the other is compressive, a transition from the phenomenological criterion to the fracture criterion is neces- 
sary. The details concerning the transition are rather complex and will not be reviewed here. Thus the remainder 
of this section outlines the details of how Powers uses Mohr-Coulomb theory to predict reliability of monolithic 
ceramic components. Mohr-Coulomb theory defines failure when the shear stress on an unspecified failure plane 
reaches a critical value. Thus the shear stress at failure is a function of the coefficient of internal friction (p), a 
compressive stress acting normal to the failure plane, and a cr (which was defined previously). The failure 
criterion takes the form 


r| + pa - a, 


(2.46) 


Powers assumes that the envelope for the Mohr-Coulomb failure criterion coincides with the envelope from 
Batdorfs fracture criterion at a = 0, which results in T = 0 cr at this point. This represents a major drawback 
and is discussed later. When viewed in Mohr’s stress space, the Mohr-Coulomb failure envelope is linear; hence 
a slope and a point completely define the envelope. In her analysis. Powers specifies the slope as the material 
parameter p. This quantity is defined as a single valued deterministic parameter. The associated point is the 
quantity o cr from Batdorfs analysis. This quantity is the only random variable in the analysis. 

Powers defines an additional quantity a max in her analysis. When viewing a state of stress in Mohr’s 
stress space, this represents the point where the failure envelope becomes tangent to Mohr s circle. In terms of 
the principal stresses and the parameter p, 0 max takes the form 


fl + P 


2 ° l ~ a2 + p CTl + ° 2 


(2.47) 


Again, index notation is briefly suspended in order to use the notation as it appears in Powers original work. 
This quantity is referred to as a "maximum" because it defines the limit of relevance from a reliability stand- 
point; that is, no reduction in reliability occurs if the failure envelope lies above this point. 


The significance of a max is easily seen when deriving the failure function. Since plane stress conditions 
simplify the analysis, again oiily the details of a surface flaw analysis are presented. In a manner similar to 
Batdorfs approach, a critical stress and orientation of a failure plane must be defined in order to formulate the 
failure function. However, it must be emphasized that the Mohr-Coulomb theory makes no allowance for infor- 
mation concerning the physics of a crack (i.e., it is not a fracture criterion). Powers expresses the failure 
function as 
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(2.48) 
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where S cr is a ratio of the critical stress (a cr ) to the maximum stress (a^) for a given state of stress; i.e., 

c = (2.49) 


The ratio S cr facilitates the numerical integration of equation (2.48). By varying this ratio from 0 to 1. a cr is 
varied from 0 to a max . This guarantees that all relevant stress states are considered in the reliability prediction. 
Also note that 63 is similar in nature to the angle to defined in equation (2.39) (the details will be discussed 
shortly), and the other terms in equation (2.48) have been defined in previous sections. 

Several combinations of the failure envelope and Mohr’s circle must be studied. Figure 2.4 shows that the 
failure envelope may not intersect Mohr’s circle, the envelope may be tangent to the circle, or the envelope may 
intersect the circle at no more than two points. To determine the location of the intersection points and the por- 
tion of the circumference of Mohr’s circle intersected. Powers formulates equation (2.46) in terms of the princi- 
pal stresses, the ratio S cr , and cos 2 (J). Here <J> locates the intersection points on the circle. The equation now 
becomes quadratic in terms of cos 2 <j>, and thus has two values, <j) j and <j> 2 . The quantity <o, which appears 
in equation (2.48), is functionally dependent on the two angles tj^ and <j> 2 in the following manner: 

_ 2(<|> l - <J> 2 ) (2.50) 

co = _1 L 

TC 


The angles c{> ^ and <J> 2 vary depending on the stress state. When there is no contact between the failure enve- 
lope and Mohr’s circle (see fig. 2.4(a)), ({Tj - $ 2 = 0, and 63 = 0^ When the failure envelope and Mohr’s circle 
are tangent at a point (fig. 2.4(b)), <t>, = <J> 2 * 0, and once again 63 = 0. When the failure envelope intercepts 
Mohr’s circle twice (fig. 2.4(c)), 63 represents the portion of Mohr’s circle that lies outside of the failure enve- 
lope. With this interpretation of 63, and the definition of ip given by equation (2.48), Powers’ model yields the 
following expression 


“k B C,£“C'‘ 1 S„dA] ( ; 

for component reliability when both principal stresses are negative. 

The Mohr-Coulomb theory is a two-parameter failure theory that represents a straight line in the Mohr 
circle stress space. The slope and a point of intersection on the t axis are used to define the linear envelope. 
Broad assumptions are made that p is deterministic and a cr is a random variable. The effects of these 
assumptions are evident when comparing compression reliability predictions to tensile reliability predictions. 
Since p is not treated as a random variable, the predicted scatter in strength (quantified by the Weibull 
modulus) is the same for simple compression and tensile tests using Powers’ model. This result has never been 
supported by experimental data for a given material. In fact, intuition says that the scatter should be greatly 
reduced in compression. Thus the Weibull modulus for compression data should be substantially higher than the 
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Weibull modulus for tensile data. Accommodating this type of behavior is not a problem with the reliability 
model presented in chapter IV. 

The following chapter presents the details of a three-parameter phenomenological failure model (the 
Mohr-Coulomb failure model is a two-parameter model). The method of transforming this deterministic phe- 
nomenological failure criterion into a reliability model is discussed in chapter IV. Basically, each of the three 
model parameters are treated as random variables with separate Weibull distributions. This leads to much greater 
flexibility in modeling reliability for multiaxial states of stress. 



Figure 2.1 . — Stress traction on an arbitrary plane. 
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Figure 2.2. — Unit sphere in three-dimensional space. 



Figure 2.3. — Intersection of Mohr's circle and Batdorf's failure 
envelope. 
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(a) No intersection between failure 
envelope and Mohr’s circle. 



(b) Envelope tangent to circle at a point. 



(c) Failure envelope cuts through Mohr's circle. 

Figure 2.4. — Power's failure envelope and 
Mohr's circle. 
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CHAPTER m 


WILLAM-WARNKE FAILURE CRITERION 

A failure criterion defines a limit state. Within this limit state, a structural component will perform its 
design application in some acceptable fashion. A primary function of the design engineer is to define what is 
acceptable performance. Performance standards depend on the design variables used to define the limit state. 
Design variables, which may include strength parameters, cyclic load limits, and allowable deformation, can be 
assembled in an n-dimensional vector 


»<, ■ (y„ Y 2 Y «) 

and the limit state function, which stipulates how the design variables interact, is expressed in general as 


" ° 


(3.2) 


This function defines a surface in the n-dimensional design variable state. If a design point (i.e., an operational 
state where each design variable has a specified value) lies within the surface, then the design point represents a 
successful operational state. If the design point falls on the surface, the component fails. For deterministic 
analyses, points outside of the failure surface are inaccessible, since failure results once the surface is reached. 

In this report, the design variable space, defined by the vector <3/ a , is limited to strength parameters for 
ceramic material systems. Since strength parameters are associated with components of the Cauchy stress tensor, 
the general functional dependence of a limit state function is expressed as 



(3.3) 


where Cfy represents the Cauchy stress tensor. A three-parameter strength criterion developed by Wiliam and 
Warnke (1975) will serve as the limit state function of primary interest here. The Willam-Warnke failure 
criterion (developed for isotropic materials) is a unified failure criterion in the sense that one limit state function 
defines failure for all regions of the stress space. As a comparison. Powers’ model adopted Batdorfs theory for 
tensile regions of the stress space and Mohr-Coulomb theory for compressive regions of the stress space. The 
Willam-Warnke criterion uses stress invariants to define the functional dependence on the Cauchy stress 
specifically 

8 K> Ip h- h) ■ 0 (3 ' 4) 


This guarantees that the function is form invariant under all proper orthogonal transformations. Here I| is the 
first invariant of the Cauchy stress J 2 is the second invariant of the deviatoric stress S y, and J 3 is the 
third invariant of the deviatoric stress. These quantities are defined as 


" °ij J ^ij°kk 


(3.5) 
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(3.6) 


h - 4 V/ 


J 3 - y 


(3.7) 


(3.8) 


where a is the identity tensor. Admitting to the function allows for a dependence on hydrostatic stress. 

The invariant J 3 allows different behavior in tension and compression since this invariant changes sign when 
the direction of a stress component is reversed. 

Wiliam and Warnke defined the limit state function with the following expression 


g 





where 



and 



(3.9) 


(3.10) 


(3.11) 


The functions B and X will be defined momentarily. The strength parameters that comprise the design vector 
ay include the uniaxial tensile strength of the material Y ( the equal biaxial compressive strength Y^, and the 
uniaxial compressive strength Y c . This model is referred to as a three-parameter model, since three strength 
parameters, - (Y t , Y c , Y^), are used to define the limit state function. Failure occurs when g - 0 and the 
multiaxial criterion is completely defined in all regions of the six-dimensional stress space. 

Since the limit state function is dependent on the six components of the Cauchy stress tensor, a design 
point and its relative position to the failure surface can be depicted in various stress spaces. Graphical 
representations can take place in a two- or three-dimensional stress space, using the components of the Cauchy 
stress tensor as coordinate axes. However, the function and the physical implications associated with the func- 
tion can be viewed completely in the three-dimensional stress space where the principal stresses serve as ortho- 
gonal coordinate axes (see fig. 3.1(a)). This space is known as the Haigh-Westergard stress space. In this 
coordinate system, a given stress state, that is, a design point PfCj, <J 2 , o 3 ), can be readily decomposed into 
hydrostatic and deviatoric components. This decomposition is shown in figure 3.1(b). Line d in figure 3.1(b) 
represents the hydrostatic axis where Oj = a 2 ” ct 3 . Point P in this stress space represents an arbitrary state of 
stress. The vector NP represents the deviatoric component of the arbitrary stress state, and the vector ON 
represents the hydrostatic component. 
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The plane passing through the origin normal to the hydrostatic line is called the n-plane. For an isotropic 
material, a failure surface projected onto the tt- plane must exhibit a sixfold symmetry. In the most general case 
where the isotropic material possesses equal strengths in tension and compression, the failure surface in the 
7 T-plane can be represented by two limiting cases (see fig. 3.2). The first case is represented by a circular failure 
surface, and the second is represented by a polygon inscribed within the circular failure surface. Any failure sur- 
face that does not fall within these two surfaces permits nonconvex regions to exist along the failure surface. 
However, proof of convexity also implies that level surfaces of a function are closed surfaces. An open region 
of the failure surface allows the existence of a load path along which failure will never occur. Thus for a con- 
vex surface, load paths cannot be traversed towards open regions of the failure surface, since open regions will 
not exist. 

A failure surface projected onto the 7 T-plane can be described conveniently with polar coordinates (r, 0). 
Here 0 is defined as an angle measured clockwise from the a r axis, and r(0) is the distance from the hydro- 
static axis to the failure surface (see fig. 3.2). Note that r(0) is a function of 0 for the inscribed polygon, and 
a constant for the circular failure surface. Physically, r(0) represents the deviatoric component of a stress state, 
since this vector lies in the 7 T-plane. In figure 3.2, the distance from the hydrostatic axis to the failure surface 
along a compressive principal axis r c is equal to the distance along a tensile principal axis q for both limiting 
cases. However, ceramic material systems exhibit very different strengths in tension and compression. Failure 
models must account for this behavior, and this can be done simply by constructing the function r(6) such that 
the intercepts along the tensile and compressive principal axes are different. The Willam-Warnke criterion 
accounts for this type of behavior by taking r c > r t (see fig. 3.3). 

As mentioned previously, isotropic materials must exhibit a sixfold symmetry in the 7T-plane. Wiliam and 
Warnke postulated that a single sector (0 < 0 < n/3) of the failure surface in the rr-plane could be represented 
as a segment of an ellipse. The major and minor axes of the ellipse were formulated as functions of the inter- 
cepts r c and r t (see fig. 3.4). Note that the minor axis of the ellipse is assumed to coincide with a tensile axis. 
However, the center of the ellipse does not necessarily coincide with the hydrostatic axis. The intercepts and 
r c depend on the strength parameters Y t , Y c , and Ybc- Equations are given later that detail the interrelation- 
ships. In general the distance r(0) is defined as 


2r c( r c 2 " r t 2 ) cos 6 + r c( 2r t “ r c) 

4c 2 - f2 ) COs2 6 + 5 ? - Vc_ 

1/2 

4 ( r c ■ r , 2 ) cos 2 6 + ( r c ■ 2r t ) 2 



where 0 <0 ti/ 3. A detailed derivation of this expression can be found in Chen (1982). Note that equa- 

tion (3.12) yields r(0) = r t for the special case of 0 = 0. Similarly r(0) = r c for 0 = 71/3. 

With the definition of r(0) given in equation (3.12), X from equation (3.9) can be expressed as 

X = (2/5)^[l/r(e)] (3 

Here X is implicitly dependent on J 3 through the angle 0. The dependence of 0 on J 3 results from the 
similarities between the trigonometric identity 
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(3.14) 


COS 3 0 - — cos 0 - _L COS ( 30 ) = 0 
4 4 


and the cubic equation used to find the deviatoric invariants 

S 3 - J 2 S - J 3 = 0 (3 ’ 

where the roots (S,, S 2 , S 3 ) are the eigenvalues of the deviatoric stress matrix. Substituting S = Y cos 0 into 
equation (3.15) yields 


, J 7 cos 0 
cos 3 0 - — 

Y 2 



By comparing this expression with equation (3.14), it becomes apparent that 

2\[h~ 

Y “ 


(3.16) 


(3.17) 


and 

cos (3 0) = ^ 3 (3.18) 

The angle 0 was first defined by Lode (1926), and the relationship between this angle and the deviatoric invar- 
iants was given by Nayak and Zienkiewicz (1972). 

Details of the derivation for the Willam-Warnke criterion have been discussed in the context of the 
7 T-plane. Since the criterion is represented by a conic surface in the three-dimensional Haigh-Westergard stress 
space, and the function is sensitive to the hydrostatic component of the stress state, details obtained from the 
two-dimensional tt- plane are not sufficient to completely describe the criterion. To gain a complete view of the 
criterion, a cutting plane is passed through the conic surface such that the entire length of the hydrostatic axis is 
contained in the cutting plane. This plane will intersect the surface along two lines. By definition, these lines are 
termed meridians. 

Meridians define the profile of the conic failure surface in the Haigh-Westergard stress space. The relative 
position of each meridian is defined by the angle 0. For the tensile meridian 0 = 0, and for the compressive 
meridian 0 = 7t/3. In the ir-plane, a compressive meridian is represented by point Q in figure 3.3, and a tensile 
meridian is represented by point T. For the Willam-Warnke criterion the meridians are linear, which is evident 
from the I, - Jj^~ stress space in figure 3.5 and equation (3.9). Since the meridians are linear, two points on 
a meridian will define its position. For the tensile meridian the two points used to determine its position are 
defined by a uniaxial tensile load path, an equal biaxial compressive load path, and their intersection with the 
meridian. For both load paths equation (3.18) yields a value of 6 = 0. Considering a uniaxial tensile load case, 
failure results when CTj (Cj = a) reaches Y t (with a 2 = a 3 = 0). The load path for this case is defined by the ratio 
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I 


l 



1.73 


The value of I, (= a) and (= o/fh) at failure will fix the position of one point on the meridian. In a 

similar fashion, the ratio 


h 



3.46 


(3.20) 


defines the load path for equal biaxial compression where ctj = a 2 (= cr) and a 3 = 0. This load path fixes the 

position of a second point on the meridian. Here Ij *= 2a and ^J 2 = o \j{2/3 ) . Failure results when a f = a 2 
= Y^. Both load paths are shown in figure 3.5. To clearly illustrate the load paths, the figure is not drawn to 
scale. The meridian is then defined by the line connecting these two points on the failure surface. 

Similarly two points are used to determine the position of the compressive meridian. The points are 
defined by a uniaxial compressive load path and the intersection of the tensile meridian with the Ij-axis. The 
load path for the uniaxial compressive case is defined by the ratio 


I, 

1 = 1.73 (3.21) 


where I,= a and = (l/fi) c. This path is shown in figure 3.5, and equation (3.18) yields a value of 
0 = n/3. Failure results when Oj = Y c with a 2 * a 3 = 0. The second point on the compression meridian is the 
tip of the failure cone. Since the tensile meridian is completely defined by the parameters Y f and Y^, the 
intersection of this line with the Ij axis provides a second point for the location of the compressive meridian. 
This point is shown in figure 3.5 as point V. The distance p from point V to the origin represents the hydro- 
static tensile stress at failure. Physically, this stress state is not easily produced in an experiment. However, this 
parameter is used to define B in equation (3.9). The parameter B is related to p by the simple expression 

B = I p (3-22) 
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As noted previously, the Willam-Warnke criterion is a three-parameter model. The parameters p, r c , and 
r, may be used to define the criterion in lieu of the strength parameters Y c> Y bc> and Y r The relationships 
are 
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Here the expressions 


and 


(6Y 1 

r - ■ w 


- v bc )’t 


~ y t 


r c = 


f 6 l 

1/2 

>’bc >’t 

l 5 J 


^bc^t + ^bc )'\ m 


(3.25) 


>’bc 


• be 


(3.26) 
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are used to simplify equations (3.23) to (3.25). Once again, the derivations of these expressions were given by 
Chen (1982). 

To gain further insight regarding physical implications of the criterion, consider the failure envelope pro- 
jected onto the a, - a 2 stress plane, which is depicted in figure 3.6. Again this is a cutting plane that passes 
through the conic surface in the Haigh-Westergard stress space. Note that, in this figure, the function defines a 
smooth failure surface for any combination of the principal stresses. Also, the differences between the tensile 
strength and compressive strength of a material are readily apparent. The ratio of the intercepts along the tensile 
and the compressive axes is equal to the ratio of Y, to Y c . This stress space is encountered again in the next 
chapter, where reliability concepts are described. 

Finally, the Willam-Warnke failure criterion degenerates to simpler models under special conditions. For 
the case of r = r. = r , where r is the same for any angle 0, the surface degenerates to a circle in the 
7 r-plane and to a cone In the three-dimensional Haigh-Westergard stress space. This is the Drucker-Prager failure 
criterion, which is a two-parameter formulation. For the special case where r c = r, = r 0 and p = °°, the model 
reduces to the single-parameter Von Mises criterion. For this case, the failure surface again becomes a circle in 
the jr-plane, but a right circular cylinder in the three-dimensional Haigh-Westergard stress space. Since " u - 
criterion exhibits no dependence on hydrostatic stress, its meridians never intersect the Ij-axis in the 7j 
stress space. 



24 



(a) Generalized failure surface in principal stress space. 



Figure 3.1. — Principal stress space. 
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Figure 3.4.— A 60° sector of Willam-Warnke failure surface depicted 
in n-plane. 



Figure 3.5. — Tensile and compressive meridians viewed in I-|_VJ 2 
stress space. 
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Figure 3.6.— Willam-Warnke failure surface depicted in o--|-a2 
stress space. 
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CHAPTER IV 


INTERACTIVE RELIABILITY MODEL 

In this chapter the details of transforming the Willam-Warnke failure function into an interactive proba- 
bilistic model are presented. Rather than predicting a fail/no fail design condition, this stochastic approach pre- 
dicts the probability of failure of a component. The design issues discussed here are similar to those outlined in 
chapter II. The interactive reliability model accounts for a reduction in reliability due to compressive stresses, 
and also accounts for decreased scatter in failure for compressive stress states in comparison to tensile stress 
states. Note that the strength parameters used in the Willam-Warnke failure criterion are treated as random vari- 
ables. Other quantities such as stiffness and loads can be treated in a probabilistic fashion (see Cruse et al., 
1988), but since the strength of ceramic-based material systems commonly varies by 100 percent or more, only 
the strength parameters are treated as random variables. 

Reliability is calculated under the assumption that the three strength parameters (Y c , Y,, and Y^) are 
independent random variables. It is assumed that each parameter is characterized by a two-parameter Weibull 
distribution (eqs. (2.1) to (2.4)); however, other distributions can be used with this approach. Using separate 
probability density functions for each random variable is versatile since other statistical distributions such as a 
three-parameter Weibull distribution or a log-normal distribution can be used to characterize the random vari- 
ables. The selection of the distribution is always dictated by the failure data. However, for the purpose of 
simplicity and illustration, only the two-parameter Weibull distribution is considered here. To define the prob- 
ability density distributions for each strength parameter, a Weibull modulus a and a scale parameter (3 must 
be determined experimentally. In general, a significant number of failure tests (a quantity which is dependent on 
the precision required for the parameter estimates of a and (3) are necessary to characterize the probability 
density function for each random variable. See Pai and Gyekenyesi (1988) for methods of parameter estimation. 

Here the functional dependence of the failure function g(^, o^) is given by equation (3.9). In general the 
reliability of a unit volume is computed from the expression 

& = Probability \g[& a , <r y ) <0 

where St represents the reliability. It is assumed that the element is homogeneous in stress; that is, no stress 
gradients exist throughout the element. Initially, the reliability calculations are based on unit volumes. Later, 
adjustments that account for arbitrary volumes are introduced. To calculate the reliability for an element of unit 
volume, the joint density must be integrated over the design space defined by the failure function. This integra- 
tion takes the form 


a - J J J Q (y c > >v ybcWc d >\ d ybc (4.2) 

where H (y c , y. , y^) is the joint density function of the random variables that correspond to the material 
strength parameters, and d s is the design space. By definition, the design space is that portion of the 
stress space bounded by the failure surface. Under the assumption that the random variables appearing in 
equation (4.2) are statistically independent, the reliability expression takes the form 
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(4.3) 


st = J J J P,(yc)P2(yt)P3(ybc) d yc d y t d ybc 


where pj(y c ), p 2 (y t ), and p 3 (y bc ) are the marginal probability density functions. These density functions are 
characterized by the two-parameter Weibull distribution. 

The integration defined by equation (4.3) yields the reliability of a unit volume. This type of integration, 
and the technique for defining the limits of integration were outlined in Sun and Yamada (1978), and 
Wetherhold (1983). To illustrate the approach, one of the simpler models presented in these references (i.e., the 
Tsai-Wu criterion) is used as an example since a closed-form solution can be obtained. Application of this type 
of criterion has been proposed for the analysis of laminate composites, since each ply is analyzed by assuming 
that the ply is an orthotropic plate subject to plane stress conditions. Here the failure function takes the form 
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where o, and o 2 are in the in-plane normal stresses, and a 6 is the in-plane shear stress. The strengths Y,, 
Y- and Y 6 are the random variables associated with the strength in the primary material direction (Y,), trans- 
verse to the primary material direction (Y 2 ), and an in-plane shear strength (Y 6 ). The reliability calculation 
follows the format outlined in equations (4.1) to (4.3), specifically. 


r = r rr 6 ’ Y2) 


Here p(y,), p(y 2 ), and p(y 6 ) are the marginal probability density functions for the respective strength variables, 
and the limits of integration are defined by systematically solving equation (4.4) for each random variable, and 
then suppressing the random variable. Thus the first limit of integration is defined by the expression 
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(4.6) 


Next, terms containing the random variable Yj are suppressed in equation (4.4), and the remainder of the 
equation is solved for Y 2 resulting in 



(4.7) 


To solve for the limits of y 6 , terms containing both Y, and Y 2 are suppressed in equation (4.4), hence 


w = o 6 


(4.8) 
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Note that each term in this particular failure function (eq. 4.4) contains one strength variable. However, for the 
Willam-Warnke model each term in the interactive formulation contains all three strength parameters. Deriving 
the limits would require substitution of equations (3.12), (3.23), (3.24), and (3.25) into equation (3.9), and 
solving for each of the strength parameters explicitly to obtain expressions similar to equations (4.6) to (4.8). 
Developing closed-form expressions is intractable because of the definition of r(0). For this reason, the triple 
integral was evaluated using Monte Carlo methods. 

The Monte Carlo technique involves generating a uniform random sample of size K for each random vari- 
able. A value is selected for each strength parameter via a random number generator. This random number is 
used with the assumed marginal probability density function (i.e., two-parameter Weibull, three-parameter 
Weibull, log normal, etc.) to obtain values for the random strength variables. Details of this computational pro- 
cedure are outlined in Wetherhold (1983). For a given stress state, the failure function is evaluated for each 
sample of random variables. Initially an element of unit volume subject to a homogeneous stress state is con- 
sidered. If g&Q. o,j) < 0 for a given trial, then that trial is recorded as a success. By repeating this process a 
suitable number of times for a given state of stress, the reliability (or cumulative distribution) of the element is 
generated. In essence, the Monte Carlo method provides a means of simulating failure experiments on a com- 
puter. This assumes that the marginal probability density functions have been suitably characterized; that is, the 
values of the a’s and the P’s are known a priori. For a sufficiently large simulation sample size, reliability is 
computed by the simple expression 



(4.9) 


where n is the number of successful trials (i.e., the number of trials where g(^ a , o^) < 0). 

Figure 4.1 shows a comparison of the Monte Carlo calculations with the underlying Weibull distribution 
assumed for the tensile strength random variable. This is a plot of probability of failure versus failure strength. 
For this case, the reliability of an element of unit volume is given by 
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where a is the applied tensile stress. For graphical purposes the natural logarithm of both sides of the 
expression is taken twice. By introducing a constant C defined as 
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then the form of equation (4.10) is 
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(4.12) 


Here the Weibull shape parameter for tensile strength defines the slope of the line depicted in figure 4.1. For 
this illustration a = 5 and P = 0.2. The three points represent estimates using the Monte Carlo method for the 
uniaxial failure strength where the specified reliabilities are 5, 50, and 95 percent. A computer algorithm which 
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numerically evaluates equation (4.5) using the Monte Carlo technique was used. For this example the desired re- 
liability, the Weibull parameters defining the marginal probability density functions for the strength parameters, 
and the desired number of Monte Carlo simulations are specified. The simulated tensile load path is described in 
chapter III. The stress state is increased incrementally along the specified load path until the desired reliabil- 
ity is found to within some predetermined error bound. The data points in figure 4.1 were generated using 100 
Monte Carlo trials. As a comparison, figures 4.2 to 4.4 show Monte Carlo estimates for different sample sizes. 
Figure 4.2 depicts predictions for 500 Monte Carlo trials; figure 4.3, 1000 trials; and figure 4.4, 10 000 trials. 
Note that as the number of trials increases, the points converge to the line representing the underlying 
Weibull distribution for the parameter Y,. This indicates that the numerical approach for evaluating equation (4.5) 
asymptotically converges to the underlying distribution as the sample size K increases. 

To illustrate behavior along other load paths, simulations were conducted for the strength parameters Y c 
and Y,^. Figure 4.5 shows the relations for the strength parameter Y c , and figure 4.6 shows the relations for 
the parameter Y^. The Weibull parameters were arbitrarily stipulated. For the uniaxial compressive case, 
a = 35 and = 2.0, and for the biaxial compressive case, a = 35 and (5 = 2.32. Note that the a’s were 
assigned higher values for the compressive strength variables than for the tensile strength variable. Although 
strength data for isotropic whisker-toughened ceramics are not available, there are sufficient experimental data 
for monolithic ceramics to indicate that compressive failure modes generally do not exhibit as much scatter as 
tensile failure modes It is believed that similar behavior will be exhibited by isotropic whisker-toughened 
materials. As in the uniaxial tensile case, estimates of reliability for 5, 50, and 95 percent were compared to the 
linear form of the Weibull equation associated with each load path. Again, all points converged with the line for 
10 000 trials. 


This technique was also used in calculating the reliability contours shown in figure 4.7. This figure 
represents the o x ~o 2 stress space. The reliability contours represent a homogeneously stressed material ele- 
ment of unit volume. Here the Weibull parameters associated with the tensile strength random variable are arbi- 
trarily chosen to coincide with the example cited in the preceding paragraph, specifically a, = 5 and ( 3 , = 0.2. 
Similarly, the Weibull parameters associated with the compressive strength random variable are arbitrarily 
specified, with a c = 35 and (5 C = 2. Finally the Weibull parameters associated with the equal biaxial com- 
pressive strength random variable are a^ = 35 and = 2.32. The three surfaces depicted in figure 4.7 
correspond to 5, 50, and 95 percent. Note that the reliability contours retain the general behavior of the 
deterministic failure surface. In general, as the a’s increase, the spacing between contours diminishes. 
Eventually with increasing a the contours would not be distinct and they would effectively map out a deter- 
ministic failure surface. An increase in the P’s shifts the relative position of the contours in an outward 
direction indicating an increase in strength. Also note that the a’s for the tensile and compressive load paths 
are different and can be specified independently of each other. This is a distinct advantage relative to the 
Powers model discussed in chapter II. Since only a tensile Weibull modulus is specified in the Powers model, 
the same scatter would occur for both tensile and compressive load paths. 

The details for computing the reliability of a single element have been presented assuming a homogeneous 
state of stress and a unit volume. To design a structural component with a varying stress field, the component is 
discretized and the stress field is characterized using finite element methods. Since component failure may 
initiate in any of the discrete elements (which typically do not have unit volumes or areas), it is useful to con- 
sider a component from a systems viewpoint. A discretized component is a series system if it fails when one of 
the discrete elements fails. This approach gives rise to weakest-link reliability theories. In a parallel system, 
failure of a single element does not necessarily cause the component to fail, since the remaining elements may 
sustain the load through redistribution. Parallel systems lead to what have been referred to in the literature as 
bundle theories. Since it is assumed that qualitatively the failure behavior of whisker-toughened ceramics 
mimics monolithic ceramics, a weakest-link reliability theory is adopted for designing structural components. 
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If the failure of an individual element is considered a statistical event, and these events are assumed to be 
independent, then the probability of failure of a discretized component is given as 

p, - 1 - n r, 13 

i=l 

where N is the number of discrete finite elements for a given component, and R, is the reliability of the I th 
discrete element. This reliability is computed in the following manner. Recall that # (the reliability based on a 
unit volume) is defined by equations (4.1) to (4.3), but calculated using the Monte Carlo techniques described 
previously. These same techniques can be used to compute Rj if the Weibull scale parameters are adjusted to 
reflect the size of the element. In general each scale parameter (P t , (3 C , and P^) is adjusted by using the 
following transformation 



(4.14) 


Here Vj is the volume of the i* element and |3* is the adjusted scale parameter. No adjustment is necessary 
for the Weibull moduli. The preceding discussion on the reliability model implied that failure of whisker- 
toughened CMC originates from volume flaws. It is quite possible that component failure is caused by surface 
and/or volume flaws; that is, competing failure modes may exist. These competing failure modes usually have 
distinctly different Weibull parameters that characterize the marginal probability density functions. Accordingly, 
equation (4.14) can be used for surface flaw analyses if Vj is replaced by the area of the ; th element, A,. 
However, for brevity, only volume flaw analysis is considered here. 

This numerical procedure has been incorporated in a public domain test bed computer algorithm given the 
acronym TCARES (Toughened Ceramics Analysis and Reliability Evaluation of Structures). Currently this 
algorithm is coupled to the MSC/NASTRAN finite element code. For a complete description of the TCARES 
algorithm, see Duffy, Manderscheid, and Palko (1989). Before using TCARES, the proposed interactive 
reliability model that has been implemented into TCARES must be characterized using an extensive data base 
that includes multiaxial experiments. It is not sufficient to simply characterize the Weibull parameters for each 
random strength variable. Multiaxial experiments should be conducted to assess the accuracy of the interactive 
modeling approach. However, once the Weibull parameters have been characterized for each random strength 
variable, the algorithm allows a design engineer to predict the reliability of a structural component subject to 
quasi-static multiaxial loads. Isothermal conditions are considered for the application that follows. However, the 
algorithm is capable of nonisothermal analyses if the Weibull parameters are specified at a sufficient and appro- 
priate number of temperature values. To illustrate certain aspects of the interactive model and the TCARES 
algorithm, a reliability analysis is performed on a test specimen which is known as the Brazilian disk. 

The Brazilian disk is used to circumvent the alignment difficulties encountered in tensile testing brittle 
materials. In addition the Brazilian disk has been used to determine tensile strengths of brittle materials that 
exhibit reduced tensile strengths relative to the compressive strength (e.g., concrete and rock). The analytical 
details concerning the stress field of the Brazilian disk have been discussed by a number of authors including 
Hondros (1959), Vardar and Finnie (1975), Chong, Smith, and Borgman (1982), and Fessler and Fricker (1984). 
Most researchers assume that tensile failure usually occurs along the diameter directly beneath the applied load 
(see fig. 4.8), splitting the disk. However, the region of the disk directly beneath the load experiences very large 
compressive stress states that dissipate slowly. The interactive model presented here allows for a reduction in 
reliability when compressive stress states are present. Thus the Brazilian disk can be used to compare the 
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interactive model with other widely used reliability models that do not account for compressive stress states. For 
simplicity, the interactive model is compared with the Principle of Independent Action (PIA) reliability model. 

It is assumed that the disk is fabricated from an isotropic whisker-toughened CMC material with a Young’s 
modulus of 300 GPa and a Poisson ratio of 0.2. A compressive pressure load of 1000 MPa was applied to the 
disk, and the subtended angle n for this example is 0.039 rad. The Weibull parameters associated with 
each random strength variable were arbitrarily chosen. Specifically the Weibull parameters associated with the 
tensile strength random variable are a, = 15 and (3, = 250. The Weibull parameters associated with the com- 
pressive strength random variable are a c = 35 and (3 C = 2500. Similarly the Weibull parameters associated with 
the equal biaxial strength random variable are a* = 35 and (3* - 2900. Note that the 0 parameters have 
units of MPa<mm) 3/a . The disk has a radius of 50 cm and a thickness of 5.0 cm, and was modeled using 1/8 
symmetry with 1044 finite elements (see fig. 4.9). The elements used in the structural analysis were 8-node 
brick elements (MSC/NASTRAN HEX/8). The tensile stress in the x-direction near the center of the disk was 
24.8 MPa. This stress remains fairly constant along the vertical diameter, except in the near vicinity of the load, 
where this stress component changes sign and becomes compressive (see fig. 4.10). The elements near the loads 
experience large compressive stresses (-997 MPa) in the y-direction that dissipate slowly down the diameter 
(see fig. 4.1 1). The stresses in this direction are compressive throughout the disk. 

When this particular discretized component was analyzed using the PIA model (with a = 1 5 and 
(3 = 250) the component reliability was 99.9 percent. Note that compressive stress states (specifically compres- 
sive principal stresses) do not affect component reliability when using the PIA model. This assumption is 
similarly adopted for other reliability theories such as Batdorf s reliability model. This lack of accounting for 
compressive stress states may be a nonconservative assumption depending on the values of the Weibull para- 
meter that characterize the compressive strength random variables. Analyzing the disk using the interactive 
reliability model presented here resulted in a component reliability of 77.7 percent. Again the Weibull para- 
meters for the compressive strength random variables were arbitrarily specified. However, a values of 35 begin 
to approach values for metals which have deterministic strength parameters, and an increase in the (3 values of 
over an order of magnitude relative to (3 t represent conservative estimates of these Weibull parameters. Thus in 
comparing the component reliability from the PIA model and the interactive reliability model, it is evident that 
accounting for compressive stress states may play an important role in the analysis of structural components. 
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Fracture stress 


Figure 4.1 . — Reliability estimates of uniaxial tensile strengths using 
100 Monte Carlo simulations. 



Fracture stress 

Figure 4.2. — Reliability estimates of uniaxial tensile strengths using 
500 Monte Carlo simulations. 
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Fracture stress 


Figure 4.3.— Reliability estimates of uniaxial tensile strengths using 
1000 Monte Carlo simulations. 



Fracture stress 


Figure 4.4. — Reliability estimates of uniaxial tensile strengths using 
1 0 000 Monte Carlo simulations. 
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Fracture stress 

Figure 4.5.— Reliability estimates of uniaxial compressive strengths 
using 10 000 Monte Carlo simulations. 



Fracture stress 

Figure 4.6.— Reliability estimates of biaxial compressive strengths 
using 1 0 000 Monte Carlo simulations. 
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Figure 4.7.— Reliability contours using the reliability model 
projected onto o-i-o^ stress space. 
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Figure 4.8. — Geometry and pressure load for Brazilian disk 
test specimen. 
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Figure 4.11. — Stresses acting on disk in normal y-direction. 
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CHAPTER V 


CONCLUDING REMARKS AND FUTURE DIRECTION 

The basic features of conducting a reliability analysis by deriving an interactive reliability model have been 
illustrated. The deterministic Willam-Warnke failure criterion serves as the theoretical foundation on which the 
reliability model was constructed. Fundamental to the work presented here is the assumption that the strength 
parameters associated with a deterministic failure criterion can be treated as random variables. As a result, the 
proposed reliability model retains the phenomenological behavior that was present in the deterministic failure 
criterion, such as sensitivity to hydrostatic stress and reduced tensile strength. The predictive capabilities 
of the interactive model were examined assuming that the two-parameter Weibull distributions characterized the 
marginal probability density functions for each random strength variable. This included both uniaxial and multi- 
axial load paths. The interactive reliability model was implemented into TCARES, a test-bed software program. 
Since this algorithm has been coupled with a general-purpose finite element program, design engineers are now 
able to use the code as a postprocessor in order to predict the reliability of a structural component subject to 
quasi-static multiaxial load conditions. A simple structural problem was presented to illustrate the reliability 
model and the computer algorithm. 

In addition, this type of reliability model can be extended to account for material anisotropy. Using ortho- 
tropic materials as an example, the parent deterministic failure function must reflect the stress state (as was done 
in this report) and the appropriate material symmetry. For orthotropy this requires that 

g = g (y, Ojj, a,, b.) (5 ' 1J 


where a ( and b ( are orthogonal unit vectors that represent the local orthotropic material directions. Because g 
is a scalar function, it must remain form invariant under arbitrary proper orthogonal transformations. Work by 
Reiner (1945), Rivlin and Smith (1969), Spencer (1971) and others demonstrated that by applying the Cayley- 
Hamilton theorem and the elementary properties of tensors, a finite set of invariants (known as an integrity 
basis) can be derived for scalar functions that are dependent on first- and second-order tensor quantities. See 
Duffy (1987) for the details regarding the application of the Cayley-Hamilton theorem for this purpose. Form 
invariance of the scalar functions is ensured if the functions depend on invariants that constitute either the 
integrity basis, or any subset thereof. A number of authors (Lance and Robinson (1971), Boehler and Sawczuk 
(1977), Arnold (1989), and Robinson and Duffy (1990)) have used this methodology to develop scalar valued 
functions that are dependent on stress (a second-order tensor) and material directions (usually characterized by 
first-order tensors as in eq. (5.1)). Clearly, the future direction alluded to here (i.e., incorporating material 
symmetry using direction tensors) is not without precedent. However, for anisotropic whisker-toughened ceramic 
composites the failure function must not only reflect the material anisotropy, but also account for reduced tensile 
strength, and a dependence on the hydrostatic component of stress, if this behavior is exhibited experimentally. 

Recall that the proposed model calculates reliability using the Monte Carlo method. For each stress state 
10 000 trials are used to compute reliability. Since this approach is used in conjunction with finite element 
methods, it could easily challenge the computational capacity of even a supercomputer as the number of discrete 
finite elements increases. In order to optimize computational efficiency, future work will also concentrate on 
using numerical schemes referred to as Fast Probability Integration (FPI) techniques. Wu (1984) outlined several 
fast probability integrators, including the methods of Rackwitz and Fiessler (1978) and Chen and Lind (1982). 
These are first-order methods since it is assumed that the limit state is linear at the design point. Quadratic 
methods have been proposed, but the added complexity is not justified by dramatic increases in accuracy. In his 
work, Wu proposed an improvement to the Rackwitz-Fiessler method. This method uses a least squares 
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technique to fit an approximated cumulative distribution function for each random variable to the true 
cumulative distribution. This approach increases accuracy with a minimal increase in computational efforts. 
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